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■(^^'bstract 
O 

present several domain decomposition algorithms for se- 
^quejitial and parallel minimization of functionals formed by 
rdiscrepancy term with respect to data and total variation 
rCnngtraints. The convergence properties of the algorithms are 
ctna'lyzed. We provide several numerical experiments, showing 
^nrfe successful application of the algorithms for the restora- 
^ tinp ID and 2D signals in interpolation/inpainting problems 
.^^pectively, and in a compressed sensing problem, for re- 
he^ering piecewise constant medical-type images from partial 
Tdujier ensembles. 



Introduction 



concrete applications for image processing, one might be 
{Jiterested to recover at best a digital image provided only 
"^fertial linear or nonlinear measurements, possibly corrupted 
noise. Given the observation that natural and man-made 
^^ages are characterized by a relatively small number of edges 
^nd extensive relatively uniform parts, one may want to help 
reconstruction by imposing that the interesting solution 
Tthe one which matches the given data and has also a few 
^continuities localized on sets of lower dimension, 
^n the context of compressed sensing [21 [8] , it has been clar- 
• ifed the fundamental role of minimizing -norms in order to 
^^^mote sparse solutions. This understanding furnishes an 
^portant interpretation of total variation minimization, i.e., 
■ the minimization of the ^i-norm of derivatives [13J, as a regu- 
larization technique for image restoration. Several numerical 
strategies to perform efficiently total variation minimization 
have been proposed in the literature. We list a few of the 
relevant ones, ordered by their chronological appearance: 

(i) the approach of Chambolle and Lions ^ by re- weighted 
least squares, see also [6 for generalizations and refinements 
in the context of compressed sensing; 
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(ii) variational approximation via local quadratic function- 
als as in the work of Vese et al. [14| [1]; 

(iii) iterative thresholding algorithms based on projections 
onto convex sets as in the work of Chambolle [3j as well as in 
the work of Combettes-Wajs [5] and Daubechies et al. [7J; 

(iv) iterative minimization of the Bregman distance as in 
the work of Osher et al. [12]; 

(v) the approach proposed by Nesterov [llj and its modifi- 
cations by Weiss et al. p!5]. 

These approaches differ significantly, and it seems that the 
ones collected in the groups iv) and v) do show presently the 
best performances in practice. However, none of the men- 
tioned methods is able to address in real-time, or at least in 
an acceptable computational time, extremely large problems, 
such as 4D imaging (spatial plus temporal dimensions) from 
functional magnetic-resonance in nuclear medical imaging, as- 
tronomical imaging or global terrestrial seismic tomography. 
For such large scale simulations we need to address methods 
which allow us to reduce the problem to a finite sequence of 
sub-problems of more manageable size, perhaps by one of the 
methods listed above. With this aim we introduced subspace 
correction and domain decomposition methods both for ii- 
norm and total variation minimizations [H [10] . Due to the 
nonadditivity of the total variation with respect to a domain 
decomposition (the total variation of a function on the whole 
domain equals the sum of the total variations on the sub- 
domains plus the size of the jumps at the interfaces), one 
encounters additional difficulties in showing convergence of 
such decomposition strategies to global minimizers. 

In this paper we review concisely both nonoverlapping and 
overlapping domain decomposition methods for total varia- 
tion minimization and we provide their properties of conver- 
gence to global minimizers. Moreover, we show numerical 
applications in classical problems of signal and image pro- 
cessing, such as signal interpolation and image inpainting. 
We further include applications in the context of compressed 
sensing for recovering piecewise constant medical- type images 
from partial Fourier ensembles [2]. 



2 Notations and preliminaries 

Since we are interested to a discrete setting, we define the 
domain of our multivariate digital signal ft = {xl < . . . < 
xj^^} X ... X {xf < . . . < xf^^} C M^, d e N and we con- 
sider the signal space H = RNixN2x...xNd ^ where Ni e N 
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for i = 1, . . . , For u e H we write u = u{x[)i^x with in- 
dex set I := Yli^iil, ' " ,Nk} and u{x[) = u{xj^ , . . . , xfj 
where ik G {l,...,A^/e} and Xi G ft. Then, for g G M^, 

the ip-noim is given by = l^fcl^) , 1 < 

j9 < oo, and for u G H the discrete gradient Vix is given 
by {\/u){xi) = ((Vii)^(xi),...,(Vii)^(xi)) with (Vii)^'(xi) = 
u{x^^ , . . . , x^^.^]^, . . . , x^^) — u{x^^ •)••••) ^i- •)••••) ^i^) if < 
and (Vii)-^(xi) = if ij = AT^-, for ah j = and for 

all i = (zi, . . . , Zfi) G X. The total variation of G in the 
discrete setting is then defined as \D{u) \ {Q) = ^ |(Vii)(xi)|, 

with 1^1 = y^ylTTTTTyj for every y = {yi,...,yd) G W^. 
We define the scalar product of ix, G as usual, (ix, v)^^ = 
and the scalar product oi p^q G as 

{p.q)h^ = EiGX^^(^i)<?n^i) + ••• Further 
we introduce a discrete divergence div : H defined, by 

analogy with the continuous setting, by div = — V* (V* is the 
adjoint of the gradient V). In the following we denote with 
ttk the orthogonal projection onto a closed convex set K. 

2.1 Projections onto convex sets 

With these notation, we define the closed convex set 



iterate 



:= u^^^ G 7Y, for example u^^^ =0, and 



(n+l) • rr/ I (^)\ 

u\ ^ argmm^^eVi J[vi ^uy) 

(n+l) _ . ^/ (n+l) , X 

^ argmm^2Gy2 J{H + ^2) 

3 Nonoverlapping domain decompo- 
sition methods 

Let us consider the disjoint domain decomposition 1] = l^i U 
1^2 and n = and the corresponding spaces Vj = {u e 
H : supp(^i) C Qj}, for j = 1, 2. Note that H = Vi ® V2. It is 
useful to us to introduce an auxiliary functional J'-^^ called the 
surrogate functional of J7: For j G {1,2} and j G {1, 2} \ {j}, 
assume a^Uj G Vj and Uj G Vj- and define 

Jl{uj+uj,a) ■.= J{uj+u.) + \\u,-af^-\\T{uj-a)f^. (3) 

As it will be clarified later, the minimization of Jj{uj a) 
with respect to Uj and for fixed u-- , a is an operation which 
can be realized more easily than the direct minimization of 
the parent functional J{uj + Uj) for the sole Uj fixed. 



K := {divp '.peW^, \pi\ < 1 for ah i e 1} . 

We briefiy recall here an algorithm proposed by Chambolle 
in [3] in order to compute the projection onto aK. The fol- 
lowing semi-implicit gradient descent algorithm is given to 
approximate 'KaK{g)'- Choose r > 0, let p^^^ = and, for any 
n > 0, iterate 



PI 



(n+l) 



pI""^ +r(V(divp^^^ 
1 + r |(V(divp(^) - g/a))i\ 



(1) 



For r > sufficiently small, the iteration adivp^^^ converges 
to TTaKig) as n — > oo (compare ^ Theorem 3.1]). 

2.2 Setting of the problem 

Given a model linear operator T : H ^ M^, we are consider- 
ing the following discrete minimization problem 

argmin{j(ii) := \\Tu - g\\l ^ 2a\Du\{n)} (2) 

where g G is a given datum and a > is a fixed regu- 
lar izat ion parameter. Note that, up to rescaling the param- 
eter a and the datum ^, we can always assume ||T|| < 1. 
Moreover, in order to ensure existence of solutions, we as- 
sume 1 ^ ker(T). For both nonoverlapping and overlap- 
ping domain decompositions, we will consider a linear sum 
H = Vi -\- V2 with respect to two subspaces Vi, V2 defined 
by a suitable decomposition of the physical domain We 
restrict our discussion to two subspaces, but the analysis can 
be extended in a straightforward way to multiple subspaces. 
With this splitting we want to minimize J' by suitable in- 
stances of the following alternating algorithm: Pick an initial 



3.1 Sequential algorithm 

In the following we denote Uj = ttvjU the orthogonal projec- 
tion onto Vj, for j = 1, 2. Let us explicitely express the algo- 
rithm as follows: Pick an initial Vi F2 3 i^i^'^^ 



,(0,M) 



G 7Y, for example u'^^^ = 0, and iterate 



y("+i.o) ^ ^{n,L) and for = 0, . . . , L - 1 

(n+1,^+1) . /Ts/ , (n,M) (n+l,. 

u\ ' = argmm^^eVi ^1 (^1 + ^ ,^1 



fn+1,0) 



(n,M) 



U2 ' = U2 ' ' and for m = 0, . . . , M — 1 
u\ ' ' = argmm^^GVs Jii^i 



^2 

,(n+l) . 



U 



(n+l,L) 



- U. 



(n+l,M) 



(n+l,m) \ 



(4) 



Note that we do prescribe a finite number L and M of inner 
iterations for each subspace respectively. 

3.2 Parallel algorithm 

The parallel version of the previous algorithm reads as follows: 



Pick an initial Vi V2 3 '^1 
example u'^^^ — 0, and iterate 



(0,L) 



,(0,M) 



/^O) G n, for 



u 



(n+1,0) _ ^ fn,L) 



U 



(n+l,£+l) 



and for £ = 0, . . . , L — 1 



argminniGVi Ji{ui 



u 



{n,M) (n+l,. 



U. 

.(n+l) . 



(n+1,0) 



(n,M) 



2 — 1^2 "'' for m 

(n+l,m+l) . (n,L) 

'2 = argminn2Gy2 ^2 (^1 



0,...,M- 1 



(n+l,m)N 



(5) 

Note that is the average of the current iteration and 

the previous one as it is the case for successive overrelaxation 
methods (SOR) in classical numerical linear algebra. 
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4 Overlapping domain decomposi- 
tion methods 

Let us consider the overlapping domain decomposition Q = 
Qi U Q2 and f^i D 7^ and the corresponding spaces Vj = 
{u G H : supp(ii) C r^j}, for j = 1,2. Note that now = 
Vi + V2 is not anymore a direct sum of Vi and F2, but just 
a Unear sum of subspaces. We define the internal boundaries 
r^- = dQj n j e {l, 2} and J e {1, 2} \ {j} (see Figured]). 

. , 

r2 

I e e 1 



— Reconstruction 



— Reconstruction 



(a) 



(b) 



^^1 

Figure 1 : Overlapping domain decomposition and internal bound- 



4.1 Sequential algorithm 

Associated to the decomposition {Qj : j = 1,2} let us fix a 
partition of unity {xj : j = 1, 2}, i.e., X1+X2 = 1, supp(xj) C 
Qj^ and Xj\rj = 0. Let us explicitely express the algorithm 



Figure 2: Here we present two numerical experiments related to 
the interpolation of a ID signal by total variation minimization, 
provided only information only out of an interval (indicated in 
green color in the figures). On the left we show an application 
of algorithm (|6]) when no correction with the partition of unity is 
provided. In this case, the sequence of the local iterations u[^'' , u^^-* 
is unbounded. On the right we show an application of algorithm 
(|6]) with the use of the partition of unity which enforces also the 
uniform boundedness of the local iterations 



now as follows: Pick an initial Vi -\-V2 3 u'l^ + : = 
7Y, for example u'^^^ =0, and iterate 



,.(0) 



of the nonoverlapping decomposition there is no additional 
linear constraint, whereas for the overlapping case we ask for 
the trace condition TjUj = Uj\rj = 0.) Such iteration can be 
explicitely computed 



u 



u 



(n+1,0) _ 
1 — 
(n+l,£+l) 



,(n) 



and for ^ = 0, . . . , L — 1 



argmin ^^e^i ^71^(1^1 

ni|ri=0 



- W 



(n) (n+1,^ 



U 



(n+l,£+l) 



fn+1,0) _ ~(n) 



nv,T%g-Tu 

(n+l 



and for m = 0, 



(n+l,m+l) 



(n+l) 



U 



u 



u. 



argmin ^^^v^ Jii^i 

n2|r2=0 
(n+l,L) , (n+l,M) 
— ""l -\-U2 

(n+l) 



,M- 1 

(n+l,L) , 



Tu 



(n+] 



(n+1,^ 
U2-, U2 



X2 ' U 



(n+l) 



(6) 

A few technical tricks are additionally required for the 
boundedness of the iterations in algorithm (j6|) with respect to 
the nonoverlapping version. First of all the local minimiza- 
tions are restricted to functions which vanish on the internal 
boundaries. Moreover since ii^^^ is formed as a sum of local 



components u{ 



(n) ,,(n) 



which are not uniquely determined on 



the overlapping part, we introduced a suitable correction by 
means of the partition of unity {xj : j = 1,2} in order to 
enforce the uniform boundedness of the sequences of the lo- 
cal iterations u^2^. With similar minor modifications, we 
can analogously formulate a parallel version of this algorithm 
as in ([5]). 

4.2 The solution of the local iterations 

The inner iterations 



)for a suitable Lagrange multiplier 7^(^+i'^) which has the role 
of enforcing the linear constraints Uj G Vj and TjUj = 0; 
^(n+i/) approximated by an iterative algorithm, see 

[TOl Proposition 4.6] for details. Note that we have to im- 
plement repeatedly the projection tTq^x for which the Cham- 
bolle's algorithm ([1]) is used. More efficient algorithms can 
also be used such as iterative Bregman distance methods [12] 
or Nesterov's algorithm [ 1 1] . 



5 Convergence properties 

These algorithms share common convergence properties, 

which are listed in the following theorem. 

Theorem. (Convergence properties) The algorithms 

and produce a sequence (u^'^^)nen with the following 

properties: 

(i)J{u^''^) > for all n e N (unless u^""^ = 

^(n+l) J. 



(a) limn^cx) 11^ 



(n+l) 



i(")||2 = 0; 



(n+l,£+l) 



argmm ^.^vj Jj [uj + u\ \ u) 



where VjUj = is some linear constraint, are crucial for the 
concrete realizability of the algorithm. (Note that in the case 



(Hi) the sequence has subsequences which con- 

verge in H; if {u''^''^)keN a converging subsequence, and 
u^"^^ is its limit, then u^"^^ is always a minimizer of J in the 
case of algorithm (|6]) (overlapping case), whereas for algo- 
rithms ([5]) (sequential and parallel nonoverlapping cases) 
this can be ensured under certain sufficient technical condi- 
tions, see Theorem 5.1 and Theorem 6.1] for details. 
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6 Numerical experiments 

In the Figure [21 Figure [3l and Figure [H we illustrate the re- 
sults of several numerical experiments, showing the successful 
application of algorithms (|4j) and (|6]), for the restoration of 
ID and 2D signals in interpolation/inpainting problems re- 
spectively, and for a compressed sensing problem. 
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Figure 3: This figure shows an application of algorithm (|6| for an 
inpainting problem. In this simulation the problem was split via 
decomposition into four overlapping subdomains. 
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Figure 4: We show an application of algorithm (|4| in a classi- 
cal compressed sensing problem for recovering piecewise constant 
medical-type images from given partial Fourier data. In this simu- 
lation the problem was split via decomposition into four nonover- 
lapping subdomains. On the top-left figure, we show the sampling 
data of the image in the Fourier domain. On the top-right we 
show the back-projection provided by the sampled frequency data 
together with the highlighted partition of the physical domain into 
four subdomains. The bottom figures present intermediate itera- 
tions of the algorithm. 
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